home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / gnu / amiga / libm.lzh / libm / README < prev    next >
Text File  |  1992-02-10  |  14KB  |  280 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  *
  33.  * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan.
  34.  * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85.
  35.  *
  36.  *    @(#)README    5.4 (Berkeley) 10/9/90
  37.  */
  38.  
  39. ******************************************************************************
  40. *  This is a description of the upgraded elementary functions (listed in 1). *
  41. *  Bessel functions (j0, j1, jn, y0, y1, yn), floor, and fabs passed over    *
  42. *  from 4.2BSD without change except perhaps for the way floating point      *
  43. *  exception is signaled on a VAX.  Three lines that contain "errno" in erf.c*
  44. *  (error functions erf, erfc) have been deleted to prevent overriding the   *
  45. *  system "errno".                                                           *
  46. ******************************************************************************
  47.  
  48. 0. Total number of files: 40
  49.  
  50.         IEEE/Makefile   VAX/Makefile    VAX/support.s   erf.c       lgamma.c
  51.         IEEE/atan2.c    VAX/argred.s    VAX/tan.s       exp.c       log.c
  52.         IEEE/cabs.c     VAX/atan2.s     acosh.c         exp__E.c    log10.c
  53.         IEEE/cbrt.c     VAX/cabs.s      asincos.c       expm1.c     log1p.c
  54.         IEEE/support.c  VAX/cbrt.s      asinh.c         floor.c     log__L.c
  55.         IEEE/trig.c     VAX/infnan.s    atan.c          j0.c        pow.c
  56.         Makefile        VAX/sincos.s    atanh.c         j1.c        sinh.c
  57.         README          VAX/sqrt.s      cosh.c          jn.c        tanh.c
  58.  
  59. 1. Functions implemented :
  60.     (A). Standard elementary functions (total 22) :
  61.         acos(x)                 ...in file  asincos.c 
  62.         asin(x)                 ...in file  asincos.c
  63.         atan(x)                 ...in file  atan.c
  64.         atan2(x,y)              ...in files IEEE/atan2.c, VAX/atan2.s
  65.         sin(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  66.         cos(x)                  ...in files IEEE/trig.c,  VAX/sincos.s
  67.         tan(x)                  ...in files IEEE/trig.c,  VAX/tan.s
  68.         cabs(x,y)               ...in files IEEE/cabs.c,  VAX/cabs.s
  69.         hypot(x,y)              ...in files IEEE/cabs.c,  VAX/cabs.s
  70.         cbrt(x)                 ...in files IEEE/cbrt.c,  VAX/cbrt.s
  71.         exp(x)                  ...in file  exp.c
  72.         expm1(x):=exp(x)-1      ...in file  expm1.c
  73.         log(x)                  ...in file  log.c
  74.         log10(x)                ...in file  log10.c
  75.         log1p(x):=log(1+x)      ...in file  log1p.c
  76.         pow(x,y)                ...in file  pow.c
  77.         sinh(x)                 ...in file  sinh.c
  78.         cosh(x)                 ...in file  cosh.c
  79.         tanh(x)                 ...in file  tanh.c
  80.         asinh(x)                ...in file  asinh.c
  81.         acosh(x)                ...in file  acosh.c
  82.         atanh(x)                ...in file  atanh.c
  83.                         
  84.     (B). Kernel functions :
  85.         exp__E(x,c) ...in file exp__E.c, used by expm1/exp/pow/cosh
  86.         log__L(s)   ...in file log__L.c, used by log1p/log/pow
  87.         libm$argred ...in file VAX/argred.s, used by VAX version of sin/cos/tan
  88.  
  89.     (C). System supported functions :
  90.         sqrt()      ...in files IEEE/support.c, VAX/sqrt.s
  91.         drem()      ...in files IEEE/support.c, VAX/support.s
  92.         finite()    ...in files IEEE/support.c, VAX/support.s
  93.         logb()      ...in files IEEE/support.c, VAX/support.s
  94.         scalb()     ...in files IEEE/support.c, VAX/support.s
  95.         copysign()  ...in files IEEE/support.c, VAX/support.s
  96.         rint()      ...in file  floor.c
  97.  
  98.  
  99.    Notes: 
  100.        i. The codes in files ending with ".s" are written in VAX assembly 
  101.           language. They are intended for VAX computers.
  102.  
  103.           Files that end with ".c" are written in C. They are intended
  104.           for either a VAX or a machine that conforms to the IEEE 
  105.           standard 754 for double precision floating-point arithmetic.
  106.  
  107.       ii. On other than VAX or IEEE machines, run the original math 
  108.           library, formerly "/usr/lib/libm.a", now "/usr/lib/libom.a", if
  109.       nothing better is available.
  110.  
  111.      iii. The trigonometric functions sin/cos/tan/atan2 in files "VAX/sincos.s",
  112.           "VAX/tan.s" and "VAX/atan2.s" are different from those in
  113.           "IEEE/trig.c" and "IEEE/atan2.c".  The VAX assembler code uses the
  114.           true value of pi to perform argument reduction, while the C code uses
  115.           a machine value of PI (see "IEEE/trig.c").
  116.  
  117.  
  118. 2. A computer system that conforms to IEEE standard 754 should provide 
  119.                 sqrt(x),
  120.                 drem(x,p), (double precision remainder function)
  121.                 copysign(x,y),
  122.                 finite(x),
  123.                 scalb(x,N),
  124.                 logb(x) and
  125.                 rint(x).
  126.    These functions are either required or recommended by the standard.
  127.    For convenience, a (slow) C implementation of these functions is 
  128.    provided in the file "IEEE/support.c".
  129.  
  130.    Warning: The functions in IEEE/support.c are somewhat machine dependent.
  131.    Some modifications may be necessary to run them on a different machine.
  132.    Currently, if compiled with a suitable flag, "IEEE/support.c" will work
  133.    on a National 32000, a Zilog 8000, a VAX, and a SUN (cf. the "Makefile"
  134.    in this directory). Invoke the C compiler thus:
  135.  
  136.         cc -c -DVAX IEEE/support.c              ... on a VAX, D-format
  137.         cc -c -DNATIONAL IEEE/support.c         ... on a National 32000
  138.         cc -c  IEEE/support.c                   ... on other IEEE machines,
  139.                                                     we hope.
  140.  
  141.    Notes: 
  142.       1. Faster versions of "drem" and "sqrt" for IEEE double precision
  143.          (coded in C but intended for assembly language) are given at the
  144.          end of "IEEE/support.c" but commented out since they require certain
  145.          machine-dependent functions.
  146.  
  147.       2. A fast VAX assembler version of the system supported functions
  148.          copysign(), logb(), scalb(), finite(), and drem() appears in file 
  149.          "VAX/support.s".  A fast VAX assembler version of sqrt() is in
  150.          file "VAX/sqrt.s".
  151.  
  152. 3. Two formats are supported by all the standard elementary functions: 
  153.    the VAX D-format (56-bit precision), and the IEEE double format 
  154.    (53-bit precision).  The cbrt() in "IEEE/cbrt.c" is for IEEE machines 
  155.    only. The functions in files that end with ".s" are for VAX computers 
  156.    only. The functions in files that end with ".c" (except "IEEE/cbrt.c")
  157.    are for VAX and IEEE machines. To use the VAX D-format, compile the code 
  158.    with -DVAX; to use IEEE double format on various IEEE machines, see 
  159.    "Makefile" in this directory). 
  160.  
  161.     Example:
  162.         cc -c -DVAX sin.c               ... for VAX D-format
  163.  
  164.        Warning: The values of floating-point constants used in the code are
  165.                 given in both hexadecimal and decimal.  The hexadecimal values
  166.                 are the intended ones. The decimal values may be used provided 
  167.                 that the compiler converts from decimal to binary accurately
  168.                 enough to produce the hexadecimal values shown. If the
  169.                 conversion is inaccurate, then one must know the exact machine 
  170.                 representation of the constants and alter the assembly
  171.                 language output from the compiler, or play tricks like
  172.                 the following in a C program.
  173.  
  174.                         Example: to store the floating-point constant 
  175.  
  176.                              p1= 2^-6 * .F83ABE67E1066A (Hexadecimal)
  177.  
  178.                         on a VAX in C, we use two longwords to store its 
  179.                         machine value and define p1 to be the double constant 
  180.                         at the location of these two longwords:
  181.  
  182.                         static long  p1x[] = { 0x3abe3d78, 0x066a67e1};
  183.                         #define      p1      (*(double*)p1x)
  184.  
  185.     Note:  On a VAX, some functions have two codes. For example, cabs() has
  186.        one implementation in "IEEE/cabs.c", and another in "VAX/cabs.s". 
  187.            In this case, the assembly language version is preferred.
  188.  
  189.  
  190. 4. Accuracy. 
  191.  
  192.             The errors in expm1(), log1p(), exp(), log(), cabs(), hypot()
  193.             and cbrt() are below 1 ULP (Unit in the Last Place).
  194.  
  195.             The error in pow(x,y) grows with the size of y. Nevertheless,
  196.             for integers x and y, pow(x,y) returns the correct integer value 
  197.             on all tested machines (VAX, SUN, NATIONAL, ZILOG), provided that 
  198.             x to the power of y is representable exactly.
  199.  
  200.             cosh, sinh, acosh, asinh, tanh, atanh and log10 have errors below
  201.             about 3 ULPs. 
  202.  
  203.             For trigonometric and inverse trigonometric functions: 
  204.  
  205.                 Let [trig(x)] denote the value actually computed for trig(x),
  206.  
  207.                 1) Those codes using the machine's value PI (true pi rounded):
  208.                    (source codes: IEEE/{trig.c,atan2.c}, asincos.c and atan.c)
  209.  
  210.                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below 
  211.                    1 ULP compared with sin(x*pi/PI), cos(x*pi/PI), and 
  212.                    atan(x)*PI/pi respectively, where PI is the machine's
  213.                    value of pi rounded. [tan(x)] returns tan(x*pi/PI) within
  214.                    about 2 ULPs; [acos(x)], [asin(x)], and [atan2(y,x)] 
  215.                    return acos(x)*PI/pi, asin(x)*PI/pi, and atan2(y,x)*PI/pi
  216.                    respectively to similar accuracy.
  217.  
  218.  
  219.                 2) Those using true pi (for VAX D-format only):
  220.                    (source codes: VAX/{sincos.s,tan.s,atan2.s}, asincos.c and
  221.                    atan.c)
  222.  
  223.                    The errors in [sin(x)], [cos(x)], and [atan(x)] are below
  224.                    1 ULP. [tan(x)], [atan2(y,x)], [acos(x)], and [asin(x)] 
  225.                    have errors below about 2 ULPs. 
  226.  
  227.  
  228.             Here are the results of some test runs to find worst errors on 
  229.             the VAX :
  230.  
  231.     tan   :  2.09 ULPs          ...1,024,000 random arguments (machine PI)
  232.     sin   :  .861 ULPs          ...1,024,000 random arguments (machine PI)
  233.     cos   :  .857 ULPs          ...1,024,000 random arguments (machine PI)
  234.     (compared with tan, sin, cos of (x*pi/PI))
  235.  
  236.     acos  :  2.07 ULPs          .....200,000 random arguments (machine PI)
  237.     asin  :  2.06 ULPs          .....200,000 random arguments (machine PI)
  238.     atan2 :  1.41 ULPs          .....356,000 random arguments (machine PI)
  239.     atan  :  0.86 ULPs          ...1,536,000 random arguments (machine PI)
  240.     (compared with (PI/pi)*(atan, asin, acos, atan2 of x))
  241.  
  242.     tan   :  2.15 ULPs          ...1,024,000 random arguments (true pi)
  243.     sin   :  .814 ULPs          ...1,024,000 random arguments (true pi)
  244.     cos   :  .792 ULPs          ...1,024,000 random arguments (true pi)
  245.     acos  :  2.15 ULPs          ...1,024,000 random arguments (true pi)
  246.     asin  :  1.99 ULPs          ...1,024,000 random arguments (true pi)
  247.     atan2 :  1.48 ULPs          ...1,024,000 random arguments (true pi)
  248.     atan  :  .850 ULPs          ...1,024,000 random arguments (true pi)
  249.  
  250.     acosh :  3.30 ULPs          .....512,000 random arguments
  251.     asinh :  1.58 ULPs          .....512,000 random arguments
  252.     atanh :  1.71 ULPs          .....512,000 random arguments  
  253.     cosh  :  1.23 ULPs          .....768,000 random arguments
  254.     sinh  :  1.93 ULPs          ...1,024,000 random arguments
  255.     tanh  :  2.22 ULPs          ...1,024,000 random arguments
  256.     log10 :  1.74 ULPs          ...1,536,000 random arguments
  257.     pow   :  1.79 ULPs          .....100,000 random arguments, 0 < x, y < 20.
  258.         
  259.     exp   :  .768 ULPs          ...1,156,000 random arguments
  260.     expm1 :  .844 ULPs          ...1,166,000 random arguments
  261.     log1p :  .846 ULPs          ...1,536,000 random arguments
  262.     log   :  .826 ULPs          ...1,536,000 random arguments
  263.     cabs  :  .959 ULPs          .....500,000 random arguments
  264.     cbrt  :  .666 ULPs          ...5,120,000 random arguments
  265.  
  266.  
  267. 5. Speed.
  268.  
  269.         Some functions coded in VAX assembly language (cabs(), hypot() and
  270.     sqrt()) are significantly faster than the corresponding ones in 4.2BSD.
  271.         In general, to improve performance, all functions in "IEEE/support.c"
  272.         should be written in assembly language and, whenever possible, should
  273.     be called via short subroutine calls.
  274.  
  275.  
  276. 6. j0, j1, jn.
  277.  
  278.         The modifications to these routines were only in how an invalid
  279.         floating point operations is signaled.
  280.